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Abstract 

We present a new method to describe the kinetics of driven lattice gases with particle-particle 
interactions beyond hard-core exclusions. The method is based on the time-dependent density 
functional theory for lattice systems and allows one to set up closed evolution equations for mean 
site occupation numbers in a systematic manner. Application of the method to a totally asymmetric 
site exclusion process with nearest-neighbor interactions yields predictions for the current-density 
relation in the bulk, the phase diagram of non-equilibrium steady states and the time evolution of 
density profiles that are in good agreement with results from kinetic Monte Carlo simulations. 

PACS numbers: 05.50.+q, 05.60. Cd, 05.70.Ln 
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Driven lattice gases are an active topic in non-equilibrium statistical mechanics due 
to both their manifold applications and their importance in fundamental studies of non- 
equilibrium systems (for reviews, see [l^). Prominent examples for applications to non- 
equilibrium processes in nature are biopolymerization jg], unidirectional motions of motor 
proteins along filaments or microtubuli p], |8| , flow of molecules through vessels or ion conduc- 
tion through membrane channels j^l, llo|. incoherent transport of electrons in molecular wires 



11] , surface growth [l2| and traffic [l3|- With respect to fundamental aspects, questions 



pertaining to the theoretical description of boundary-induced phase transitions and of non- 



equi 



ibrium steady states (NESS) can be studied within a conceptually simple framework 



|3j. 

A standard model of a driven lattice gas system is the asymmetric simple exclusion 
process (ASEP), which refers to the directional stochastic hopping transport of particles on 
a one-dimensional lattice with hard-core exclusions that prevent multiple occupation of a 
lattice site. In an open system with particle injections and ejections at the left and right 
boundaries, this model exhibits phase transitions between non-equilibrium steady states of 

nn 

different mean site occupation in the bulk [14J, [15[. Analytical methods like the matrix 
approach or Bethe ansatz have been developed and successfully applied to calculate exactly 
the distribution of microstates in the NESS 2M5|. Various extensions of the ASEP have 



3een considered in the past, as, for example, particle rods covering several 



161 ] . or stochastic injection/ejection (Langmuir 



however, comparatively few studies exist 
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attice sites 



dnetics) of particles in the bulk 1^] . So far, 
19| which are concerned with the description 



of driven lattice gas systems, where the particles do not only interact via athermal hard-core 
exclusion. 

In this work we will focus on developing a method for treating the kinetics of such 
systems with interactions beyond hard-core exclusions. In general, starting with the master 
equation for the Markovian time evolution of the microstates, exact evolution equations 
can be derived for the mean site occupation numbers (henceforth called densities). These 
have the form of a continuity equation with currents depending on equal-time correlations 
of occupation numbers (henceforth called "correlators"). The challenge is to develop proper 
methods to calculate these correlators so that closed evolution equations for the densities 
result that well account for the kinetic behavior. In the treatment of interaction effects in 
the driven system considered by Schiitz et al. [19| this problem was addressed by utilizing 



the special jump rates introduced in the Katz-Lebowitz-Spohn model [181 ]. which yield a site 
occupation statistics in a driven bulk system that can be exactly mapped onto an equilibrium 
system with nearest-neighbor interactions (Ising system). As a consequence, the current- 
density relation in the bulk could be determined exactly and, by employing the maximum 
and minimum current principle Ijl, [19J], the phase diagram of NESS in an open system 
calculated in good agreement with Monte Carlo simulations. 

In this Letter we show how by using the time-dependent density functional theory 
(TDFT) for lattice systems a systematic method is provided that allows one to treat the 
kinetics of driven lattice gases with interactions. To demonstrate the procedure, we consider 
a totally asymmetric site exclusion process (TASEP) with nearest-neighbor interactions. 
Predictions of the theory for the current- density relation in the bulk, the phase diagram 
of NESS, the density profiles in the NESS, and the time evolution of density profiles are 
compared to results of kinetic Monte Carlo (KMC) simulation and yield a surprisingly good 
agreement for this far-from-equilibrium system. 

The Hamiltonian for a one-dimensional lattice gas with nearest-neighbor interactions V 

is 

H (n) = V ^ n i n i+i ■ C 1 ) 

i 

The set of occupation numbers n = {rij} specifies the microstate of the system, where n« = 
or 1 if site i is vacant or occupied by a particle, respectively (nf = rti). For unidirectional 
nearest-neighbor hopping considered in the TASEP, the stochastic dynamics of the system 
is specified by the rates Ti(n) for a particle on site i to jump to a vacant neighboring site 
(z + 1) in the configuration n. Starting from the master equation for the time evolution of the 
probability density P(n,t) of microstates, one derives the discrete version of the continuity 
equation for the densities pi{t) = (rij) t (see, e. g. [20| for a general derivation), 

^=J,-At)-Jdt), (2) 

where the mean currents from i to % + 1 are given by 

jiif) = (m (1 - n i+l ) Yi{n)) t , (3) 

and (•••)( denotes an average with respect to P(n,t). To complete the specification of the 
dynamics we define the jump rates Ti(n) as 22] 

r-i(n) = uexp [- (U(n^) - H(n)) /2k B T] . (4) 



Here v is an attempt frequency, and n^'* +1 ^ refers to the target configuration of the jump, 
where, with respect to the initial configuration n, rii and rii+i are interchanged, while the 
other rik are the same. We use v~ x as time unit and k^T as energy unit in the following 
(z/ -1 = 1 and k-^T = 1). With Eq. (jl]) the currents can be written explicitely as 

ji = {hi-xnin i+1 h i+ 2)t + e v/2 \n i ^ 1 n i hi + ihi + 2)t + e~ v/2 (hi-inin i+1 ni +2 )t + (n i _in i n i+ in i+2 >t , 

(5) 

where we introduced the hole occupation numbers hi — 1 — rii. 

In order to express the correlators in Eq. (jSJ) in terms of the densities, we now apply 
the TDFT, which is based on the (time-)local equilibrium approximation. According to 
this, P(n,t) is expressed by the Boltzmann probability associated with %{n) plus a time- 
dependent external potential U[n, p(t)] = £\ Vi[p(t)]rii that equals the potential which 
would generate p(t) = {pi{t)} as equilibrium density profile (according to the classical 
version of the Mermin theorem this potential is unique for given interaction). In effect this 
implies that the correlators at any time t are supposed to be related to the densities as in 
an equilibrium system. 

To find the equilibrium correlator-density relations, we apply the methods developed in 



the Markovian approach to derive exact density functionals 2l| . Accordingly we express the 
joint probabilities p^ (rii, . . . , n i+ j) for the occupation numbers rij, . . . , n i+ j in equilibrium 
by the Markov chain pi{ +1) (n,, . . . , n i+j ) = p^(ni)Y^ s=1 w(ni + s\n i+s -i), where piq(nj) is 
the equilibrium probability for rii, and w(rii + i\rii) = (ni,n i+ i)/p^q (rii) the conditional 
probability for rij + i if rii is given. Since the joint probabilities are directly connected to the 
correlators, e. g., piq (rii_i = 0, = 1, n i+1 = 0, n i+2 = 1) = (^i-i^i^i+i^j+2)eq, all correlators 
involving more than two occupation numbers in Eq. ([5]) can thus be reduced to two-point 
correlators. The TDFT expression for the current hence becomes 

j7 dft = [{ Pl+2 - C l+l )e- y ' 2 + ~ Pl+l - p l+2 + Ci+J - + e^CU] , (6) 

PiPi+i 

with p 



via 



2l| 



I — p,i and Ci = (niTii+i)t- The two-point correlators are related to the densities 

Ci = exp(-V) (Pi-Ci)(Pi + i-C i ) t (7) 
J- — Pi — Pi+i + 

which can be explicitly solved to yield functions Ci = Ci(pi, Pi+i). In this way, the currents 
jTDFT are gj ven as functionals of the density profile pit). 



To test the quality of the TDFT approach we start by considering a homogeneous system 
with periodic boundary conditions. In this case we suppress the site indices in Eq. (J7J) and set 
ji = j(p) in Eq. (EJ). For V — > 0, j(p) approaches the parabola j = p — p 2 for particles feeling 
only hard-core repulsion. When V exceeds a critical value V c = —2 In (v5 — 2) ~ 2.89, j(p) 
develops a double-hump structure with two maxima at densities pi 2 = Pi 2^V)- I n the hmit 
V ->■ oo we find p\ 2 = 1/2 ± (y/2 - l) /2, in agreement with earlier findings reported by 
Krug jl^j], and j — > (x 3 / 2 — 2x + x 1 / 2 ) /(1—x) with x = (2p — l) 2 , meaning that the particle 
movement is frozen for a half-filling system. 

In Fig. [JJwe compare the analytical findings (a) for the pair correlator C(p) and (b) the 
current- density relation j(p) with results from KMC simulations for various V. For V < V c 
as well as for large V > AV C , we find an excellent agreement with the simulation data, see 
also the inset in Fig. Ufa), which shows the interaction dependence of C(p — 0.5). For 
intermediate interaction strength V c < V < 4V C and near half-filling, some deviations occur. 
An interesting effect is seen in Fig. [Hb) for weak couplings V < V c : both the theoretical 
predictions and the simulations show an increased current compared to the case V = 0. This 
phenomenon is caused by effective particle-hole attraction for weak coupling strengths. 

Next we consider an open system of N sites that is coupled to two particle reservoirs L 
and R at the left and right boundaries with densities pl and respectively. The dynamics 
of injection (ejection) of particles from the left (right) reservoir is defined such that the same 
correlator-density relations apply at the boundaries as in the corresponding bulk (periodic 



ring) systems with densities p^ and p&. We used the method described in [19[ to implement 
this model in the KMC simulations. 

Let us first see how well the TDFT approach captures the time evolution of density 
profiles. To this end we consider reservoir densities pl = 0.4 and pr = 0.1 and an initially 
empty lattice system at time t — 0. For the interaction we choose deliberately V = 2V C , 
where comparatively large deviations were seen in Fig[TJ Such a "worst-case" choice permits 
the best evaluation of approximation limits. In Fig. [2] we compare the results from numerical 
solutions of the nonlinear rate equations 02]) [where the ji are given by Eqs. and (0)] with 
the KMC results for four different times and in the stationary limit t — > oo. The excellent 
agreement between the two data sets demonstrates the power of our TDFT approach for 
driven lattice gas systems with interactions. In the inset of Fig. [2] we show the steady state 
(t — > oo) profile for the corresponding non- interacting case V = 0. As one would expect 



from an application of the maximum and minimum current principle to the bulk current- 
density relation shown in Fig. QJb), different bulk densities pe are found for the two different 
interaction strengths: Pb = Pl = 0.4 for V = and ps = p\{2V c ) ~ 0.31 for V = 2V C . 

More generally, using the minimum and maximum current principle, we can evaluate the 
phase diagram of the NESS for arbitrary interaction strengths based on the bulk current- 
density relation derived above. Alternatively we can use the t — > oo limit of the rate 
equations to identify the bulk densities of the NESS. As expected, we find that the results 
of both procedures agree 23j. As an example we show in Fig. El^a) the phase diagram of 
the NESS for V = 2V C , and for each of the occurring seven phases we display the stationary 
density profile in Fig. |3]^b). In all cases the predictions are well confirmed by the KMC 
results. The largest deviations occur for the transition lines between the phases I /III, I/VII 
and the "mirror lines" IV/II and IV/VII. Note that the phase diagram in Fig. [3(a) exhibits 
particle-hole symmetry as required by the Hamiltonian in Eq. (CQ) and the dynamics specified 
in Eq. (j3J). There are two maximal current phases V and VI characterized by the bulk 
densities pe = p* and pe = p£, respectively, and one minimal current phase with pe = 0.5. 
The remaining four phases are determined by the reservoir densities pl and pr. 

In summary, we presented a new analytical approach to describe the kinetics of driven 
lattice gas systems with interactions beyond hard-core exclusions based on the TDFT. The 
approach was demonstrated for a TASEP with nearest-neighbor interaction and provided 
very good results for the kinetic behavior. It is clear that the TDFT is not restricted to 
the TASEP situation or nearest-neighbor interactions, but can be applied also to ASEPs 
(or purely boundary driven systems) as well as other types of interactions. In the TDFT 
closed nonlinear evolution equations for the densities result due to local relations between 
correlators and densities, which allow one to capture the kinetics much better than simple 
factorization schemes (note that simple factorization schemes of correlators would completely 
fail for V > V c in the example studied in this work). Although one can regard it as a 
weakness, it is, on the hand, advantageous that the TDFT relies on density functionals for 
equilibrium systems. The development and improvement of such equilibrium functionals in 
condensed matter systems have been an intensive research area in the past, in particular also 
with respect to find good approximation schemes for the setup of functionals in dimensions 
larger than one. Hence one can expect that these developments will be useful to treat more 
complex non-equilibrium systems. We hope that our findings will stimulate further research 

6 



in this direction. 

We thank W. Dieterich for very valuable discussions. 
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FIG. 1. (Color online) (a) Two-point correlator C{p) and (b) current-density relation j(p) for 
various V calculated by Eqs. (|6|) and ([7]) (lines). The inset in Fig. Uta) shows the pair correlator 
as a function of V/V c for fixed density p = 0.5. The analytical results are compared with the true 
non-equilibrium quantities obtained by Monte Carlo simulations (circles) for a ring system with 
1000 sites. Each data point results from an average over a time corresponding to 10 9 particle jumps 
in the steady state. 
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FIG. 2. (Color online) Time evolution of pi from an initially empty lattice with pl = 0.4, /jr = 0.1, 
V = 2V C and = 1000. TDFT solutions are marked by lines. The KMC results (circles) are 
averaged over 10 6 configurations. For V = 0, the steady-state density profile is shown in the inset. 
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FIG. 3. (a) Steady-state phase diagram and (b) density profiles of a system with 1000 sites and 

V = 2V C . The phase diagram contains seven phases which are represented by density profiles 
corresponding to I (p B = PL, low-density phase): pl = 0.1, pr = 0.4; II (p B = PL,)'- PL = 0.6, 
Pr = 0.7; III (p B = Pr): Pl = 0.3, p R = 0.4; IV (p B = Pr, high-density phase: p L = 0.6, p R = 0.9; 

V (pb = p*, maximal-current phase): pl = 0.6, pr = 0.1; VI (ps = p%, maximal-current phase): 
PL = 0.9, pr = 0.4; VII (pb = 0.5, minimal-current phase): pl = 0.3, pr = 0.7. Solid and dashed 
lines in the phase diagram, which denote first- and second-order phase transitions, respectively, 
as well as density profiles are calculated using TDFT. Data points are obtained by Monte Carlo 
simulations with the same statistics as in Fig. [TJ 



11 



